%% Processing of rice pile data
%First half of script processes the raw time series and plots a power spectra

% Input raw sediment feeder data 
close all
clear all
fid = fopen('sediment_scale_log_03172022.txt'); %input mass log .txt file 
out = textscan(fid,'%s %f %f','delimiter',',','headerlines',1);
fclose(fid);
datetime = out{1};
mass = out{2}; %mass data 
dt1 = datetime(1);%time experiment started
sec = etime(datevec(datetime),datevec(dt1)); %calulates a time series
maxsec = floor(sec(end));

%Interpolate raw data onto time series with regular dt
dt = 1.12;%~dt between mass measurements
t = 1:dt:maxsec; %time into run in seconds
t = t';
mass_int = interp1(sec,mass,t); %interploate mass data onto regular time series
M = movmedian(mass_int,5);
Q = (M(2:end) - M(1:end-1))/dt;%Q = mass flux in grams/sec
Q = [0;Q];
T = max(t); %max time
nt = size(M,1);%total number of mass measurements
Qm = mean(Q);
Qmax = max(Q);

figure(1)
plot(t,Q,'k')
xlabel('Time (s)')
ylabel('Sediment flux (grams/s)')
axis([0 max(t) 0 max(Q)])
taxis = max(t);

% Probability of exceedance from efflux data 
Q2 = Q(Q>0.1); %Use data from Q which is >0.1 grams/sec 
nbins = 10.^[-1:0.1:1.8]; %Define bin edges logarithmically - -1:0.1:1.8, 0.1g:1.25g:63g
[N,edges] = histcounts(Q2,nbins,'Normalization','count'); %Count how many points in Q2 fall between the defined bin ranges 
n = length(Q2);
N_s = N/n;
P = zeros(1,length(N_s)) %Calculate for each bin, the number of measurements >= to each bin of intrest
for i = 1:length(N_s)
    P(i) = sum(N_s(i:end));
end  
c = nbins(:,1:end-1) + nbins(:,2:end)/2 %Find the center of each bin

figure(2)
loglog(c,P,'ok')
xlabel('Qs (grams/sec)')
ylabel('Probability of exceedence')

% Centering of data and generation of n ergodic realizations from total dataset
Qdetrended = Q - mean(Q);%Centering of data
ne = 1;%number of ergodic realizations from full data set
nte = floor(nt/ne);
Qdetrended_e = [];
for i = 1:ne
    qdetrended_e = Qdetrended(1+i*nte-nte:i*nte);
    Qdetrended_e = [Qdetrended_e qdetrended_e];
end
Qdetrended_e = Qdetrended_e';

% Plot MTM spectra 
addpath('E:\MATLAB scripts\Experimental confidence bands\Rice_Realizations\Rice_Realizations') %add path to where confidence band file is
nw = 4;%number of windows
nsim = 100;%number of random simulations based on AR(1) model fit, to figure out confidence bands
t1 = t(1:nte);
t1 = t1';
MTMe = [];
CHI95 = [];
for i = 1:ne;
    datax = Qdetrended_e(i,:);
    datax = [t1' datax'];
    [fd1, po, theored1, red90, red95, red99,red90global,red95global, red99global] ...
   = bpl(datax, nw, 0, 0.5, 0, 1, 1, 'grams');
    MTMe = [MTMe po];
    CHI95 = [CHI95 red95];
end
MTMmean = mean(MTMe');%ensemble average data
CHImean = mean(CHI95');
if ne == 1;%just to deal with data if you aren't actually ensemble averaging anything.
    MTMmean = MTMe;
    CHI95 = CHI95;
end
period = 1./fd1;
lnperiod = log(period);%convert data over to log space
lnMTMmean = log(MTMmean);%convert data over to log space
dp = 0.001;%pick an increment that is fine, relative to data spacing, to interpolate unevenly spaced data onto.
lnperiod_int = min(lnperiod):dp:max(lnperiod(2:end));%create vector with spacing dp to interpolate onto.
lnMTMmean_int = interp1(lnperiod(2:end),lnMTMmean(2:end),lnperiod_int);%Do the interpolation
mean_power = 4465 %mean(exp(lnMTMmean_int));
MTMmean = MTMmean./mean_power;

figure(3)
loglog(period,MTMmean,'k','linewidth',1);
axis([min(period) max(period) min(MTMmean) max(MTMmean)])
ylabel('Spectral Power (Qs)')
xlabel('Period (s)')

% Quantify autogenic timescales
figure(4)
loglog(period,MTMmean,'r','linewidth',2);
ylabel('Spectral Power')
xlabel('period (sec)')
axis([min(period) max(period) min(MTMmean) max(MTMmean)])
hold on
lnperiod = log(period);
lnMTMmean = log(MTMmean);
dp = 0.001;
lnperiod_int = min(lnperiod):dp:max(lnperiod(2:end));
lnMTMmean_int = interp1(lnperiod(2:end),lnMTMmean(2:end),lnperiod_int);
FC = findchangepts(lnMTMmean_int,'Statistic','linear','MaxNumChanges',2);
C1 = [exp(lnperiod_int(FC(1)+1));exp(lnperiod_int(FC(1)+1))];
C2 = [exp(lnperiod_int(FC(2)+1));exp(lnperiod_int(FC(2)+1))];
YC = [min(MTMmean) max(MTMmean)];
loglog(C1,YC,'--k')
loglog(C2,YC,'--k')


%% Processing cyclic experiments
% Second half of script has processing for determining signal degradation
% and signal detectability

% Calculate signal degradation

% Reshape Q into sections equal to the signal period and average across the row
Qe = reshape(Q(1:29004),12,[]); % if array isnt divisible by the period round data to nearest value
Qe_mean = mean(Qe,2); %take mean across the row
time = 1:length(Qe_mean);

% input flux sine wave 
period = 12;
amp = 0.37;
mean_flux = 0.37;
ps = 0; %phase shift
wave = amp*sin((2*pi/period)*(time+ps))+mean_flux;

x1 = time'; %transpose time
y = Qe_mean;
yu = max(y);
yl = min(y);
yr = (yu-yl);                % Range of ‘y’
yz = y-yu+(yr/2);
zx = time(yz .* circshift(yz,[0 1]) <= 0); % Find zero-crossings
P = period;                   %Period of signal
ym = mean(Qe_mean);                % Estimate offset
fit = @(b,x)  b(1).*(sin(2*pi*x./P + 2*pi/b(2))) + b(3);      % Function to fit
fcn = @(b) sum((fit(b,x1) - y).^2);                              % Least-Squares cost function
s = fminsearch(fcn, [yr;  -1;  ym])                       % Minimise Least-Squares                
xp = linspace(min(x1),max(x1),P);

figure(5)
plot(time,wave,'r','Linewidth',1)
xlabel ('Time (s)')
ylabel('Mean ensemble sediment flux (grams/sec)')
hold on
plot(time,Qe_mean,'k','Linewidth',1)
hold on
plot(xp,fit(s,xp), 'b','Linewidth',1)
legend('Input signal','Ensemble mean efflux','Measured output signal')
legend('Location','eastoutside')

s = s';
amp_ratio = (s(1)/amp)*100 % amplitude difference

%Calculate signal detectability 

